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We numerically investigate the mechanical properties of static packings of ellipsoidal particles in 
2D and 3D over a range of aspect ratio and compression A<f>. While amorphous packings of spherical 
particles at jamming onset (A(j> = 0) are isostatic and possess the minimum contact number Zi so 
required for them to be collectively jammed, amorphous packings of ellipsoidal particles generally 
possess fewer contacts than expected for collective jamming (z < Zi so ) from naive counting argu- 
ments, which assume that all contacts give rise to linearly independent constraints on interparticle 
separations. To understand this behavior, we decompose the dynamical matrix M = H — S for 
static packings of ellipsoidal particles into two important components: the stiffness H and stress 
S matrices. We find that the stiffness matrix possesses N(zi so — z) eigenmodes eo with zero eigen- 
values even at finite compression, where N is the number of particles. In addition, these modes eo 
are nearly eigenvectors of the dynamical matrix with eigenvalues that scale as A(j>, and thus finite 
compression stabilizes packings of ellipsoidal particles. At jamming onset, the harmonic response of 
static packings of ellipsoidal particles vanishes, and the total potential energy scales as S 4 for per- 
turbations by amplitude S along these 'quartic' modes, eo. These findings illustrate the significant 
differences between static packings of spherical and ellipsoidal particles. 

PACS numbers: 83.80.Fg61.43.-j,63.50.Lm,62.20.-x 
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I. INTRODUCTION 

There have been many experimental [IHH], computa- 
tional [3411, an d theoretical studies of the structural 
and mechanical properties of disordered static packings 
of frictionless disks in 2D and spheres in 3D. In these 
systems, counting arguments, which assume that all par- 
ticle contacts give rise to linearly independent impenetra- 
bility constraints on the particle positions, predict that 
the minimum number of contacts required for the sys- 
tem to be collectively jammed is N c > iV™ ln = N^of + 1, 
where N^of — Nd for fixed boundary conditions and 
^Vdof = Nd — d for periodic boundary conditions P, , 
where d is the spatial dimension and N is the number 
of particles The additional contact is required be- 

cause contacts between hard particles provide only in- 
equality constraints on particle separations [9]. In the 
large-system limit, this relation for the minimum num- 
ber of contacts reduces to z > Zi so , where z = 2N C /N 
is the average contact number. Disordered packings of 
frictionless spheres are typically isostatic at jamming on- 
set with z = zjso, and possess the minimal number of 
contacts required to be collectively jammed @. Fur- 
ther, it has been shown in numerical simulations that 
collectively jammed hard-sphere packings correspond to 
mechanically stable soft-sphere packings in the limit of 
vanishing particle overlaps @, [12, EH • 

In contrast, several numerical and experimen- 

tal studies [18|, [l9| have found that disordered packings 
of ellipsoidal particles possess fewer contacts (z < z; so ) 
than predicted by naive counting arguments, which as- 



sume that all contacts give rise to linearly independent 
constraints on the interparticle separations. Despite this, 
these packings were found to be mechanically stable 

(ms) uKU. 

In a recent manuscript fl4| by Donev, et ai, the au- 
thors explained this apparent contradiction — that static 
packings of ellipsoidal particles are mechanically stable, 
yet possess z < z iso . The main points of the argument 
are included here. The set of N c interparticle contacts 



imposes N c constraints, 



1 



fij/^ij < 0, where r# 



is the center-to-center separation and Oij is the contact 
distance along fy between particles i and j. In disor- 
dered MS sphere packings with z — Zj SO , each of the N c 
interparticle contacts represents a linearly independent 
constraint. In contrast, some of the N c contacts for MS 
packings of ellipsoidal particles give rise to linearly de- 
pendent constraints. Linearly dependent constraints do 
not block the degrees of freedom that appear in the con- 
straint equations for sphere packings, whereas they can 
block multiple degrees of freedom in packings of ellip- 
soidal particles because they have convex particle shape 
with a varying radius of curvature [Til ]. 

For static packings of spherical particles, interparti- 
cle contacts give rise to only "convex" constraints (Fig.[T] 

(a) ), while contacts can yield "convex" or "concave" con- 
straints in packings of ellipsoidal particles (Figs. [T] (a) and 

(b) ). Note that the distinction between convex and con- 
cave constraints is different than the distinction between 
convex and concave particles. For example, ellipsoids 
have a convex particle shape, but static ellipsoid pack- 
ings can possess concave interparticle constraints. 

In jammed packings of ellipsoidal particles, there are 
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FIG. 1: Schematic diagram that illustrates locally (a) con- 
vex (positive radius of curvature) and (b) concave (negative 
radius of curvature) interparticle constraints in packings of 
hard ellipsoidal particles. Inaccessible regions (with fij > 0) 
are shaded blue. The axes labeled ui and U2 indicate two 
representative orthogonal directions in configuration space. 



N(zi so — z) special directions eo in configuration space 
along which perturbations give rise to interparticle over- 
laps that scale quadratically with the perturbation am- 
plitude, fij ~ S 2 0, E3- Displacements in all other di- 
rections yield overlaps that scale linearly with S, fij ~ S, 
as found for jammed sphere packings. This novel scal- 
ing behavior for packings of ellipsoidal particles can be 
explained by decomposing the dynamical matrix M — 
H — S for these packings into two important compo- 
nents: the stiffness matrix H that contains all second- 
order derivatives of the total potential energy V with 
respect to the configurational degrees of freedom, and 
the stress matrix S that includes all first-order deriva- 
tives of V with respect to the particle coordinates. The 
directions eo are the eigenvectors of the stiffness matrix 
H with zero eigenvalues. 

For static packings of ellipsoidal particles at the jam- 
ming threshold (A<j> = 0) that interact via purely repul- 
sive linear spring potentials (i.e. V ~ ffj), we find that 
the total potential energy increases quartically when the 
system is perturbed by S along the eo directions, V ex c5 4 , 
where the constant c > 0. Also, at the jamming thresh- 
old, the stress matrix S = and zero modes of the stiff- 
ness matrix are zero modes of the dynamical matrix. 

In this manuscript, we will investigate how the me- 
chanical stability of packings of ellipsoidal particles is 
modified at finite compression (A<fi > 0). For example, 
when a system at finite A<p is perturbed by amplitude 5 
along eo, do quadratic terms in S arise in the total poten- 
tial energy or do the contributions remain zero to second 
order? If quadratic terms are present, do they stabilize 
or destabilize the packings, and how do the lowest fre- 
quency modes of the dynamical matrix scale with A(f> 
and aspect ratio? 

We find a number of key results for static packings of 
ellipsoidal particles at finite compression (A<p > 0) in- 
cluding: 1) Pac kings o f ellipsoidal particles generically 
satisfy z < Zi so [l#4l7| : 2) The stiffness matrix H pos- 
sesses N(z- lso — z) eigenmodes eo with zero eigenvalues 



even at finite compression (A(f> > 0); and 3) The modes 
eo are nearly eigenvectors of the dynamical matrix (and 
the stress matrix —S) with eigenvalues that scale as cA<j>, 
with c > 0, and thus finite compression stabilizes pack- 
ings of ellipsoidal particles [15|. In contrast, for static 
packings of spherical particles, the stiffness matrix con- 
tributions to the dynamical matrix stabilize all modes 
near jamming onset. At jamming onset (A<f> — 0), the 
harmonic response of packings of ellipsoidal particles van- 
ishes, and the total potential energy scales as <5 4 for per- 
turbations by amplitude 5 along these 'quartic' modes, 
eo- Our findings illustrate the significant differences be- 
tween amorphous packings of spherical and ellipsoidal 
particles. 

The remainder of the manuscript will be organized as 
follows. In Sec. [TT1 we describe the numerical methods 
that we employed to measure interparticle overlaps, gen- 
erate static packings, and assess the mechanical stability 
of packings of ellipsoidal particles. In Sec. IIII1 we de- 
scribe results from measurements of the density of vibra- 
tional modes in the harmonic approximation, the decom- 
position of the dynamical matrix eigenvalues into contri- 
butions from the stiffness and stress matrices, and the 
relative contributions of the translational and rotational 
degrees of freedom to the vibrational modes as a func- 
tion of overcompression and aspect ratio using several 
packing-generation protocols. In Sec. IIV1 we summarize 
our conclusions and provide promising directions for fu- 
ture research. We also include two appendices. In Ap- 
pendix [XJ we show that the formation of new interparti- 
cle contacts affects the scaling behavior of the potential 
energy with the amplitude of small perturbations along 
eigenmodes of the dynamical matrix. In Appendix [Bj 
we provide analytical expressions for the elements of the 
dynamical matrix for ellipse-shaped particles in 2D that 
interact via a purely repulsive linear spring potential. 



II. METHODS 



In this section, we describe the computational meth- 
ods employed to generate static packings of convex, 
anisotropic particles, i.e. ellipses in 2D and prolate el- 
lipsoids in 3D with aspect ratio a = a/b of the major to 
minor axes (Fig. and analyze their mechanical prop- 
erties. To inhibit ordering in 2D, we studied bidisperse 
mixtures (2-to-l relative number density), where the ra- 
tio of the major (and minor) axes of the large and small 
particles is ai/a s = bi/b s = 1.4. In 3D, we focused on a 
monodisperse size distribution of prolate ellipsoids. We 
employed periodic boundaries conditions in unit square 
(2D) and cubic (3D) cells and studied systems sizes in the 
range from N = 30 to 960 particles to address finite-size 
effects. 
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FIG. 4: Ellipses with a = 2 positioned at the Gay-Berne con- 
tact distance a?j. For two ellipses with the same size, the 
(a) end-to-end configuration is exact, while the (b) side-to- 
end configuration has a 5% relative error. For two ellipses 
with dj/ai — 1.4, the (c) end-to-end configuration has a rel- 
ative error of 1%, while the (d) side-to-end configuration has 
a relative error of 10%. 




FIG. 2: We focus on (a) ellipses in 2D with aspect ratio a — 
a/b defined as the ratio of the major to minor axis and (b) 
prolate ellipsoids in 3D where a is the ratio of the polar to 
equatorial lengths. 




FIG. 3: Definition of the contact distance Oij for ellipsoidal 
particles i and j with unit vectors (a and fij that characterize 
the orientations of their major axes, cnj is the center-to-center 
separation at which ellipsoidal particles first touch when they 
are brought together along fij at fixed orientation. 



A. Contact distance 

In both 2D and 3D, we assume that particles interact 
via the following pairwise, purely repulsive linear spring 
potential 




ij < (Tij 



(1) 



scales will be expressed in units of e, I = y/I/m, and 
ly/m/e, respectively, where m and I are the mass and 
moment of inertia of the ellipsoidal particles. 

Perram and Wertheim developed an efficient method 
for calculating the exact contact distance between ellip- 
soidal particles with any aspect ratio and size distribu- 
tion in 2D and 3D [28j . In their formulation, the contact 
distance is obtained from 
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where e is the characteristic energy of the interaction, ry 
is the center-to-center separation between particles i and 
j, and oy is the orientation-dependent center-to-center thus — <Ty(A m i n ), /? = /3(A m in), X — x(A m in), and 



The approximation cr° = Cij(A = 1/2) is equivalent 
to the commonly used Gay-Berne approximation for the 
contact distance [13, [H[ . The accuracy of the Gay-Berne 
approximation depends on the relative orientation of the 
two ellipsoidal particles, and in general is more accurate 
for monodisperse systems. For example, in Fig. 01 we 
show a% for several relative orientations of both monodis- 
perse and bidisperse systems. The relative deviation from 
the true contact distance can be as large as e ~ 10% for 
dj/ai = 1-4 and a = 2. Thus, the Gay-Berne approxi- 
mation should be used with caution when studying poly- 
disperse packings of ellipsoidal particles. For monodis- 
perse ellipses with a = 2, 0% < e < 5%. We find sim- 
ilar results for 3D systems. Unless stated otherwise, we 
employ the exact expression for contact distance, and 



separation at which particles i and j come into contact 
as shown in Fig. |3J Below, energies, lengths, and time 



°ij ~ ^(Ainin); where A m i n is the minimum obtained 
from Eq. [2] 
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B. Packing generation protocol 



We employ a frequently used isotropic c omp ression 
method for soft, purely repulsive particles [30l [3lj to 
generate static packings of ellipsoidal particles at jam- 
ming onset (A0 = 0). Static packings at jamming on- 
set are characterized by infinitesimal but nonzero total 
potential energy and pressure. This isotropic compres- 
sion method consists of the following steps. We be- 
gin by randomly placing particles at low packing frac- 
tion (</> = 0.2) with random orientations and zero ve- 
locities in the simulation cell. We successively com- 
press the system by small packing fraction increments 
6<f) = 1CP 3 , with each compression followed by conjugate 
gradient (CG) energy minimization until the total poten- 
tial energy per particle drops below a small threshold, 
V/N < Vtoi = 10~ 16 , or the total potential energy per 
particle between successive iterations of the minimization 
routine is {V t +\ — V t )/V t < V to i- The algorithm switches 
from compression to decompression if the minimized en- 
ergy is greater than 2V to \- Each time the algorithm tog- 
gles from compression to decompression or vice versa, the 
packing fraction increment is halved. 

The packing-generation algorithm is terminated when 
the total potential energy per particle satisfies Vtoi < 
V/N < 2Vto\- Thus, using this method we can locate 
the jammed packing fraction cj>j and particle positions at 
jamming onset for each initial condition to within 10~ 8 . 
After jamming onset is identified, we also generate con- 
figurations at specified A</> = <j) — <pj over six orders of 
magnitude from 10 -8 to 10 -2 by applying a prescribed 
set of compressions with each followed by energy mini- 
mization. 

To determine whether the accuracy of the energy min- 
imization algorithm affects our results (see Sec. IIII Dp . 
we calculate the eigenvalues of the dynamical matrix as 
a function of the total kinetic energy (or deviation from 
zero in force and torque balance on each particle) at each 
Acj). To do this, we initialize the system with MS pack- 
ings from the above packing-generation algorithm and 
use molecular dynamics (MD) simulations with damping 
terms proportional to the translational and rotational ve- 
locities of the ellipsoidal particles to remove excess kinetic 
energy from the system [32|. The damped MD simula- 
tions are terminated when the total kinetic energy per 
particle is below K/N — Ktolj where K to \ is varied from 
10 -16 to 10 -32 . This provides accuracy in the particle 
positions of the energy minimized states in the range from 
1(T 8 to 1CT 16 . 



For the damped MD simulations, we solve Newton's 
equations of motion (using fifth-order Gear predictor- 
corrector methods 33]) for the center of mass position 
and angles that characterize the orientation of the long 



axis of the ellipsoidal particles. In 2D, we solve 

i>3 



m lm = 2^ Fi i~ brVi 



dt 2 

d 2 9i 
dt 2 



i>j 



(3) 



(4) 



where 6i is the angle the long axis of ellipse i makes 
with the horizontal axis, Hi is the translational velocity 
of particle i, 0i is the rotational speed of particle i, b r 
and bg are the damping coefficients for the position and 
angle degrees of freedom, and the moment of inertia / = 
m(a 2 + 6 2 )/4. The force F^j on ellipse i arising from an 
overlap with ellipse j is 
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where 
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-dV(r l0 )/dr l0 = (T-Ml 



-) for the purely repulsive 



linear spring potential in Eq. [TJ and f # and ipij are illus- 
trated in Fig. [5] 

To calculate the torque = \r\- x Fij] ■ z in Eq. SI 
we must identify the point of contact between particles i 
and j, 



where 



tanr. 



2 y/ a~ 2 + tan 2 T io 

[(cos 9i — sin 0i tan ry )x + 

(sin#i tan Ty + cos#j)$] , 



tan(-0ij - 0i) 
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ij 
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and r\p ipij, and are depicted in Fig. [5] From Eqs.[5] 
and[8j we find 



Ti 
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+ a~ 4 tan 2 r^) (l + a~ 2 tan 2 Ty) 
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FIG. 5: Geometry of two ellipses in contact. The angles 9i 
and 6j characterize the orientation of particles i and j relative 
to the horizontal axis, i.e. p4 — cos 0i£ + sin 8iy. tpij gives the 
angle between the center-to-center separation vector fij and 
the horizontal axis and ipij = — sin tpij -x + cos tpijy is the angle 
unit vector in polar coordinates. The unit vector Fij = —Fji 
points in the direction of the force on particle i due to particle 
j at the point of contact, ffj points from the center of particle 
i to the point of contact with particle j, and Tij is the angle 
between fa and fy. 



C. Dynamical matrix calculation 

To investigate the mechanical properties of static pack- 
ings of ellipsoidal particles, we will calculate the eigen- 
values of the dynamical matrix and the resulting de nsity 
of vibrational modes in the harmonic approximation [34] . 
The dynamical matrix is defined as 



derivative involving angles, Fg = —dV(rij)/89i, where 

2 



A 



\x[^) {2aA{B ++ B_)+ X C{Bl-Bl)), (12) 
ijij cos 9i - Xij sinfli 



a(xij cos 9i + jjij sin Oi) ± a 1 (x-i 



B± 

C = cos 2 (0i - 



Uij sin 9 j) 



(l + Xcos[^ -9j])r i3 



Complete expressions for the matrix elements of the dy- 
namical matrix for ellipses in 2D are provided in Ap- 
pendix |Bj In 3D, we calculated the first derivatives of V 
with respect to the particle coordinates analytically, and 
then evaluated the second derivatives for the dynamical 
matrix numerically. 

The vibrational frequencies in the harmonic approxi- 
mation can be obtained from the Ndf—d nontrivial eigen- 
values rrii of the dynamical matrix, u>i = y/ mi / eb s . d of 
the eigenvalues are zero due to periodic boundary condi- 
tions. For all static packings, we have verified that the 
smallest nontrivial eigenvalue satisfies m m i n /N > 10 -10 . 

Below, we will study the density of vibrational fre- 
quencies D(u) = (N(oj + Aw) - N(uj))/(N doi Auj) as a 
function of compression A(j) and aspect ratio a, where 
N((jj) is the number of vibrational frequencies less than 
uj. We will also investigate the relative contributions 
of the translational and rotational degrees of freedom 
to the nontrivial eigenvectors of the dynamical matrix, 



rrii 

for ellipses 
r 7—1 7 
i m xi 



f 7 — 1 7 — 1 7 
\ m xi ' m -" ' m 
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j= N 3- 
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:N 
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} for prolate ellipsoids in 3D, where i labels 



the eigenvector and runs from 1 to Ndf - 
eigenvectors are normalized such that rhf = 1. 



The 



Mm = 



d 2 V 
dukdui 



(11) 



where Uk (with k = 1, . . . ,dfN) represent the dfN 
degrees of freedom in the system and df is the number of 
degrees of freedom per particle. In 2D df = 3 with u — 
{xi,x 2 ,- ■ .,x N ,yi,y 2 ,. ■ -,VNhQ\h&2,- ■ -M^n} and in 3D 
for prolate ellipsoids df = 5 with u — {x\,X2,- ■ -,XN,yi, 
?/2,. ■ .,yN,zi,z 2 ,. . .,z N ,lg9i,lg92,- ■ .,lg 9N,h^>i,h4>2,. ■ ., 
h4>N}, where 0, is the polar angle and 4>i is 
the azimuthal angle in spherical coordinates, 

h = Va 2 + 6 2 /2, / 3 = .^/(a 2 + 6 2 )/5, and 

'(26 2 + (a 2 -6 2 )sin 2 

The dynamical matrix requires calculations of the first 
and second derivatives of the total potential energy V 
with respect to all positional and angular degrees of free- 
dom in the system. The first derivatives of V with respect 
to the positions of the centers of mass of the particles f\ 
can be obtained from Eq.[5] In 2D, there is only one first 



D. Dynamical matrix decomposition 

The dynamical matrix (Eq. 1111) can be decomposed 
into two component matrices M = H — S: 1) the stiff- 
ness matrix H that includes only second-order derivatives 
of the total potential energy V with respect to the con- 
figurational degrees of freedom and 2) the stress matrix 
S that includes only first-order derivatives of V. The kl 
elements of H and S are given by 



Hki — 



Ski 



E 



d 2 v d(nj j 'aij) d(nj j (Tij) 



i>3 



dinj/cTij) 2 du k 



dm 



dV 



d 2 (rij/(Tij) 



d(r l j/a l j) dukdui 



(13) 
(14) 



where the sums are over distinct pairs of overlapping par- 
ticles i and j. Since 9 2 V r /9(r i j/cr i j) 2 = e for the purely 
repulsive linear spring potential (Eq.QJ, the stiffness ma- 
trix depends only on the geometry of the packing (i.e. 
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d(rij/<Jij)/duk)- Also, at zero compression Acj) = 0, 
S = 0, M = H , and only the stiffness matrix contributes 
to the dynamical matrix. The frequencies associated with 
the eigenvalues hi of the stiffness matrix (at any Acj)) 
are denoted by com — y/hi/eb s , and the stiffness matrix 
eigenvectors are normalized such that hf = 1. 

E. Contact number 

When counting the number of interparticle contacts 
N c , we remove all rattler particles (defined as those with 
fewer than d+ 1 contacts) and do not include the contacts 
that rattler particles make with non-rattler particles [35j j . 
Removing these contacts may cause non-rattler particles 
to become rattlers, and thus this process is performed 
recursively. Note that for ellipsoidal particles with d+1 
contacts, the lines normal to the points (or planes in 3D) 
of contact must all intersect, otherwise the system is not 
mechanically stable. The number of contacts per particle 
is defined as z = N C /(N — N r ), where N r is the number 
of rattlers. We find that the number of rattler particles 
decreases with aspect ratio from approximately 5% of the 
system at a = 1 to zero for a > 1.2 in both 2D and 3D. 

III. RESULTS 

Static packings of ellipsoidal particles at jamming on- 
set typically possess fewer contacts than predicted by 
isostatic counting arguments [lij], z < Zj SO , over a wide 
range of aspect ratio as shown in Fig. [6l This finding 
raises a number of important questions. For example, 
are static packings of ellipsoidal particles mechanically 
stable at finite Acj) > 0, i.e. does the dynamical ma- 
trix for these systems possess nontrivial zero-frequency 
modes at Acj) > 0? In this section, we will show that 
packings of ellipsoidal particles are indeed mechanically 
stable (with no nontrivial zero-frequency modes) by cal- 
culating the dynamical, stress, and stiffness matrices for 
these systems as a function of compression Acj), aspect 
ratio a, and packing-generation protocol. Further, we 
will show that the density of vibrational modes for these 
systems possesses three characteristic frequency regimes 
and determine the scaling of these characteristic frequen- 
cies with Acf> and a. 

A. Density of vibrational frequencies D(uj) 

A number of studies have shown that amorphous 
sphere packings are fragile solids in the sense that the 
density of vibrational frequencies (in the harmonic ap- 
proximation) -D(w) for these systems possesses an excess 
of low-frequency modes over Debye solids near jamming 
onset, i.e. a plateau forms and extends to lower frequen- 
cies as Acj) -> | HH, H3- In this work, we will cal- 
culate D{uj) as a function of Acj) and aspect ratio a for 
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a 

FIG. 6: Average contact number z versus aspect ratio a for 
static packings of (a) bidisperse ellipses in 2D and (b) prolate 
ellipsoids in 3D at jamming onset. The isostatic values Zi so = 
6 (2D) and 10 (3D) are indicated by dashed lines. 



amorphous packings of ellipsoidal particles and show that 
the density of vibrational modes for these systems shows 
significant qualitative differences from that for spherical 
particles. 

In Figs. [7] (a) and (b), we show -D(w) on linear and 
log-log scales, respectively, for ellipse-shaped particles in 
2D at Acj> = 10 -8 over a range of aspect ratios from 
a = 1 to 2. We find several key features in D(w): 1) 
For low aspect ratios a < 1.05, D(ui) collapses with that 
for disks (a = 1) at intermediate and large frequencies 
0.25 < oj < 2.25; 2) For large aspect ratios a > 2, D(uj) is 
qualitatively different for ellipses than for disks over the 
entire frequency range; and 3) A strong peak near uj = 
and a smaller secondary peak at intermediate frequencies 
(evident on the log-log scale in Fig.[7](b)) occur in D(oj) 
for a > 1. Note that at finite compression Acj) > 0, we 
do not find any nontrivial zero-frequency modes of the 
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FIG. 7: (a) The density of vibrational frequencies D(uu) for 
N = 240 ellipse-shaped particles at A(j> = 10 -8 with aspect 
ratio a = 1.0 (solid), 1.001 (dotted), 1.05 (dashed), and 2.0 
(dot-dashed). D(uj) for a — 1 has been scaled by 2/3 relative 
to those with a > 1 to achieve collapse at low aspect ratios, 
(b) D{uj) for the same aspect ratios in (a) on a log-log scale. 
The inset illustrates the three characteristic frequencies oji, 
L02, and u>3 in D(uj) for a — 1.001. 



dynamical matrix in static packings of ellipses and ellip- 
soids. The only zero-frequency modes in these systems 
correspond to the d constant translations that arise from 
periodic boundary conditions and zero-frequency modes 
associated with 'rattler' particles with fewer than d + 1 
interparticle contacts. 

To monitor the key features of D(ui) as a function of 
A(f> and a, we define three characteristic frequencies as 
shown in the inset to Fig. [7J (b). u)% and 0J2 identify the 
locations of the small and intermediate frequency peaks 
in D(uj), and W3 marks the onset of the high-frequency 
plateau regime in D(u>). For our analysis, we define uj^ as 
the largest frequency (< 1) with D(u>) < 0.15, which is 
approximately half of the height of the plateau in -D(w) 
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FIG. 8: (a) The density of vibrational frequencies D(lj) for 
JV = 512 prolate ellipsoids at A(j> = 10~ 6 for a = 1.0 (solid), 
1.001 (dotted), 1.005 (dashed), and 1.2 (dot-dashed). D(u) 
for a — 1 has been scaled by 3/5 relative to those with a > 1 
to achieve collapse at low aspect ratios, (b) D(oS) for the same 
aspect ratios in (a) on a log-log scale. The inset illustrates 
the three characteristic frequencies wj., 0J2, and 0J3 in D(to) 
for a = 1.001. 



at large frequencies. All three characteristic frequencies 
increase with aspect ratio. Note that we only track u>2 
and W3 for aspect ratios where W2 < W3. For example, 
the intermediate and high-frequency bands characterized 
by L02 and uj^ merge for a > 1.2. 

As shown in Fig. [SJ -D(w) for 3D prolate ellipsoids dis- 
plays similar behavior to that for ellipses in 2D (Fig. [7]) 
for aspect ratios a < 1.5. For example, D(ui) for el- 
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FIG. 9: The density of vibrational frequencies D(lo) for N = 
240 ellipses as a function of compression A<j) = 10" 7 (solid), 
10" 5 (dashed), 10~ 3 (dotted), and 1CT 2 (dot-dashed) for (a) 
a = 1.05 and (b) 2. 
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FIG. 10: (a) Characteristic frequencies ui\ (circles), CJ2 
(squares), and ljs (diamonds) from D(uj) as a function of as- 
pect ratio a - 1 for N = 240 ellipses in 2D at A(j) = 10" 8 . 
The solid (dashed) lines have slope 1/2 (1). (b) ^i/(A^) 1/2 
for systems with = 240 ellipses in 2D at A<f> = 10 -7 (cir- 
cles), 10 -6 (squares), 10 -5 (diamonds), 10 -4 (upward trian- 
gles), 10 -3 , (downward triangles), and 10 -2 (crosses). The 
solid line has slope 1/2. 



lipsoids possesses low, intermediate, and high frequency 
regimes, whose characteristic frequencies u>i, 0J2, and W3 
increase with aspect ratio. Note that the intermediate 
and high-frequency bands W2 and w 3 merge for a > 1.02, 
which occurs at lower aspect ratio than the merging of 
the bands in 2D. Another significant difference is that in 
3D D(lj) extends to higher frequencies at large aspect 
ratios (a > 1.2) than D(ui) for ellipses. 

We note the qualitative similarity between the D(lo) 
for a = 1.005 ellipsoids shown in Fig.[8](b) and D(uj) for 
a = 0.96 presented in Fig. 1 (c) of Ref. [l6| for tu > 10~ 2 . 
However, Zeravcic, et al. suggest that there is no weight 
in D(ui) for uj < 10 ~ 2 except at uj = for both oblate 
and prolate ellipsoids, in contrast to our results in Fig. HI 

In Fig. IH1 we show the behavior of D(uj) for ellipse 



packings as a function of compression for two aspect 
ratios, a — 1.05 and 2. We find that the low-frequency 
band (characterized by uj{) depends on A<f>, while the in- 
termediate and high frequency bands do not. The inter- 
mediate and high frequencies bands do not change sig- 
nificantly until the low-frequency band centered at lu\ 
merges with them at A<j) s» 10~ 3 and ~ 10~ 4 for a = 1.05 
and 2, respectively. 

We plot the characteristic frequencies wi, 0J2, and 0J3 
versus aspect ratio a — 1 for ellipse packings in Fig. [10] 
and ellipsoid packings in Fig. [IT] The characteristic fre- 
quencies obey the following scaling laws over at least two 
orders of magnitude in a — 1 and five orders of magnitude 
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FIG. 11: (a) Characteristic frequencies wi (circles), 012 
(squares), and (diamonds) from D(uj) as a function of 
aspect ratio a — 1 for N = 240 prolate ellipsoids in 3D at 
A(j> = 10~ 6 . The solid (dashed) lines have slope 1/2 (1). (b) 
cji/(A</>) 1//2 for systems with N = 240 prolate ellipsoids at 
Acj> — 10~ 6 (circles), 10~ 5 (squares), and 10~ 4 (diamonds). 
The solid line has slope 1/2. 



in A<j>: 

wi ~ (A0) 1 / 2 (a - l) 1 ^ 2 , (15) 

w 2 ~ (a-1), (16) 

W3 ~ (a-l) 1/2 . (17) 

Similar results for the scaling of L02 and W3 with a — 1 
were found in Ref. jig . We will refer to the modes in 
the low-frequency band in D(lo) (with characteristic fre- 
quency u>i) as 'quartic modes', and these will be discussed 
in detail Sec. IIII CI The scaling of the quartic mode fre- 
quencies with compression, uj\ ~ (A^) 1 / 2 , has important 
consequences for the linear response behavior of ellip- 
soidal particles to applied stress [To] ]. 
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FIG. 12: The distribution of frequencies D(uuh) associated 
with the eigenvalues of the stiffness matrix H for N = 240 
ellipse packings as a function of compression A<f> = 10 -5 (dot- 
ted), 10" 3 (dashed), and 10" 2 (dot-dashed) for a = 1.05. The 
vertical solid line indicates the 'zero-frequency' tolerance Wtoli 
which is the lowest frequency obtained for the dynamical ma- 
trix for packings at a = 1.05 and the smallest compression 
(A<j>= 10" 8 ) in Fig. □ 



B. Dynamical Matrix Decomposition 

As shown in Fig. [BJ static packings of ellipsoidal parti- 
cles can possess z < Zi so over a wide range of aspect ratio, 
yet as described in Sec. IIII A[ the dynamical matrix M 
contains a complete spectrum of Nd f — d nonzero eigen- 
values ?Tij near jamming. To investigate this intriguing 
property, we first calculate the eigenvalues of the stiff- 
ness matrix iJ, show that it possesses N z 'zero'-frequency 
modes whose number matches the deviation in the con- 
tact number from the isostatic value, and then identify 
the separate contributions from the stiffness and stress 
matrices to the dynamical matrix eigenvalues. 

In Fig. 1121 we show the distribution of frequencies 
D{ujh) associated with the eigenvalues of the stiffness ma- 
trix for ellipse packings at a — 1.05 as a function of com- 
pression Acf). We find three striking features in Fig. [T2] 1) 
Many modes of the stiffness matrix exist near and below 
the zero-frequency threshold (determined by the vibra- 
tional frequencies of the dynamical matrix at a = 1.05 
and Acf> = 10 -8 ), 2) Frequencies that correspond to the 
low-frequency band characterized by uj\ are absent, and 
3) The nonzero frequency modes (with uj^ > 10~~) do 
not scale with A<p as pointed out for the dynamical ma- 
trix eigenvalues in Eqs. 1161 and 1171 Further, we find that 
the number of zero-frequency modes N z of the stiffness 
matrix matches the deviation in the number of contacts 
from the isostatic value (A^ so — N c ) for each A(j> and as- 
pect ratio. Specifically, N z — N iso — N c over the full 
range of A<j> for 99.96% of the more than 10 3 packings 
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FIG. 13: The two contributions to the dynamical matrix 
eigenvalues, (a) H and (b) — S, plotted versus cj 2 = H — S 
for ellipse packings with a = 1.05 and A<j> = 10 -6 (circles), 
10~ 5 (squares), 10~ 4 (diamonds), and 10~ 3 (triangles). In (a) 
and (b), the solid lines correspond to H — oj 2 and S — cu 2 . 
In the main panel and inset of (a), only modes correspond- 
ing the intermediate and high frequency bands are included. 
In the main panel and inset of (b), only modes correspond- 
ing the low frequency band are included. The insets to (a) 



and (b), which plot -S/(A(j>Y 



co 2 /A(j), show the deviations uP" — H = - 
intermediate- frequency modes and up" 
for low- frequency modes. 



iuj 2 and n/{A<j}) 2 



S oc A(f> for high- and 
- {-S) =Hcc (A(j>) 2 



for aspect ratio a < 1.1 and for 100% of the more than 
10 3 packings for a > 1.1. 

In Fig. [T31 we calculate the projection of the dynam- 
ical matrix eigenvectors fhi onto the stiffness and stress 
matrices, H — rh\Hrrii and S = rh\Srrii, where rh\ is the 
transpose of m, and uif = m\Mrhi = T~L — S. Fig. Q2] 
(a) shows that for large eigenvalues uij of the dynamical 
matrix (i.e. within the intermediate and high frequency 



bands characterized by W2 and 013 in Fig. [7]), the eigenval- 
ues of the stiffness and dynamical matrices are approxi- 
mately the same, H « ujf . The deviation ujf — % = — S 7 
shown in the inset to Fig.[T3](a), scales linearly with A<f>. 
Thus, we find that the intermediate and high frequency 
modes for packings of ellipsoidal particles are stabilized 
by the stiffness matrix H. 

In the main panel of Fig. [T3] (b) , we show that for fre- 
quencies in the lowest frequency band (characterized by 
uji) the eigenvalues of the stress and dynamical matrices 
are approximately the same, — S ~ ujf. In the inset to 
Fig. Q2] (b), we show that the deviation uif — (-S) = H 
scales as (A<f>) 2 . Thus, we find that the lowest frequency 
modes for packings of ellipsoidal particles are stabilized 
by the stress matrix —S over a wide range of compression 
A(f>. Similar results were found previously for packings 
of hard ellipsoidal particles [14| . In contrast, for static 
packings of spherical particles, the stress matrix contri- 
butions to the dynamical matrix are destabilizing with 
— S < for all frequencies near jamming, and H stabi- 
lizes the packings as shown in Fig. [T4l 



C. Quartic modes 



We showed in Sec. IIII Al that the dynamical matrix 
M for packings of ellipsoidal particles contains a com- 
plete spectrum of Ndf — d nonzero eigenvalues m, for 
A<fi > despite that fact that z < Zj S0 . Further, we 
showed that the modes in the lowest frequency band 
scale as u)\ ~ (A^) 1 / 2 in the A</) — ► limit. What hap- 
pens at jamming onset when A</> = 0, i.e. are these 
low-frequency modes that become true zero-frequency 
modes at A<f> = stabilized or destabilized by higher- 
order terms in the expansion of the potential energy in 
powers of the perturbation amplitude? 

To investigate this question, we apply the following 
deformation to static packings of ellipsoidal particles: 



u = uq + 5rhi , 



(18) 



where 5 is the amplitude of the perturbation, m, is an 
eigenvector of the dynamical matrix, and uq is the point 
in configuration space corresponding to the original static 
packing, followed by conjugate gradient energy minimiza- 
tion. We then measure the change in the total potential 
energy before and after the perturbation, AV. 

We plot AV/N versus 5 in Fig. [15] for perturbations 
along eigenvectors that correspond to the smallest non- 
trivial eigenvalue m\= wj^in °f ^ ne dynamical matrix for 
static packings of (a) disks and ellipses and (b) spheres 
and prolate ellipsoids at A<j> = 10 -7 . As expected, for 
disks and spheres, we find that AV/N rj mujf nin 5 2 over a 
wide range of S in response to perturbations along eigen- 
vectors that correspond to the smallest nontrivial eigen- 
value. In contrast, we find novel behavior for AV/N 
when we apply perturbations along the eigendirection 
that corresponds the lowest nonzero eigenvalue of the 
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FIG. 14: The two contributions to the dynamical matrix 
eigenvalues, (a) H and (b) ~S/ A(j>, plotted versus the eigen- 
values of the dynamical matrix to 2 for bidisperse disk packings 
at A(f> = 10 -6 (circles), 10 -5 (squares), 10 -4 (diamonds), and 
1CP 3 (triangles). In (a) the solid line corresponds to H = cu 2 . 
Note that —5 < over the entire range of frequencies. 



5 



o 




-5 -4 

log 10 6 



-3 




log 10 6 



FIG. 15: The change in the total potential energy AV (nor- 
malized by the particle number N) before and after apply- 
ing the perturbation in Eq. [18] with amplitude 8 along the 
eigenvector that corresponds to the lowest nontrivial eigen- 
value of the dynamical matrix for packings of (a) disks (solid 
line) and ellipses with a = 1.1 (dashed line) and (b) spheres 
(solid line) and prolate ellipsoids with a = 1.1 (dashed line) 
at A(f> = 10~ 7 . The dot-dashed (dotted) lines have slope 2 
(4). 



dynamical matrix for packings of ellipsoidal particles. In 
Fig. [T5] we show that AV/N obeys 



AV m 9 r9 
N 2 k 



c k S 4 



(19) 



where uik oc A^ 1 / 2 and the constants Ck > 0, for pertur- 
bations along all modes k in the lowest frequency band of 
D(u) for packings of ellipsoidal particles when we do not 
include changes in the contact network following the per- 
turbation and relaxation. (See Appendix [XI for measure- 
ments of AV/N when we include changes in the contact 
network.) Eigenmodes in the lowest frequency band are 
termed 'quartic' because at A<fi = they are stabilized 
by quartic terms in the expansion of the total potential 



energy with respect to small displacements [Fq l . 

For S Sf , the change in potential energy scales as 
AV/N - uj 2 k S 2 , whereas AV/N - c k S 4 for S > 5%, 
where the characteristic perturbation amplitude 5^ = 
Wfe\/m/2cfc. In the insets to Fig. [16] (a) and (b), we 
show that the characteristic perturbation amplitude av- 
eraged over modes in the lowest frequency band scale as 
S* ~ (A0) 1 / 2 / (a — l) 1 / 4 for static packings of ellipses in 
2D and prolate ellipsoids in 3D, which indicates that the 
Cfe possess nontrivial dependence on aspect ratio a. 

The quartic modes have additional interesting features. 
For example, quartic modes are dominated by the rota- 
tional rather than translational degrees of freedom. We 
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FIG. 16: The change in the total potential energy AV/N 
for perturbations along the 'quartic' modes (solid) and all 
other modes (dashed) as a function of amplitude S for (a) 
ellipses and (b) prolate ellipsoids with a = 1.1 for Ac/> = 
10~ 7 . The dotted (dot-dashed) lines have slope 2 (4). The 
solid vertical lines indicate the characteristic amplitude 8* at 
which AV/N crosses over from quadratic to quartic scaling 
averaged over the quartic modes. The insets show the scaling 
of S* /(Ac/>) 1//2 with a — 1 for several values of compression: 
A(j> = 10 -8 (circles) and 10 -7 (squares) for 2D systems and 
10 -6 (diamonds), 10~ 5 (triangles), and 10 -4 (pluses) for both 
2D and 3D systems. The solid lines in the insets have slope 
-0.25. 



identify the relative contributions of the translational and 
rotational degrees of freedom to the eigenvectors of the 
dynamical matrix in Figs. [T7] and 1181 The contribution 
of the translational degrees of freedom to eigenvector rhi 
is defined as 



Nd f 
3 = 1 A 



where the sum over A includes x and y in 2D and x, y, 
and z in 3D and the eigenvectors are indexed in increas- 
ing order of the corresponding eigenvalues. Since the 
eigenvectors are normalized, the rotational contribution 
to each eigenvector is Ri = 1 — T s ; . 

For both ellipses in 2D and prolate ellipsoids in 3D, we 
find that at low aspect ratios (a < 1.1), the first N (2N) 
modes in 2D (3D) are predominately rotational and the 
remaining 2N (3N) modes in 2D (3D) are predominately 
translational. In the inset to Figs. HTT b) and [TBI we show 
that T increases as (a — 1)^, where C ~ 1-5 (1-25) for el- 
lipses (prolate ellipsoids) , for both the low and intermedi- 
ate frequency modes. For a > 1.2, we find mode-mixing, 
especially at intermediate frequencies, where modes have 
finite contributions from both the rotational and trans- 
lational degrees of freedom. For a < 1.2, the modes be- 
come increasingly more translational with increasing fre- 
quency. For a > 1.2, the modes become more rotational 
in character at the highest frequencies. Our results show 
that the modes with significant rotational content at low 
a correspond to modes in the low and intermediate fre- 
quency bands of D(uS), while the modes with significant 
translational content at low a correspond to modes in 
the high frequency band of D{uj). 



D. Protocol dependence 

We performed several checks to test the robustness and 
accuracy of our calculations of the density of vibrational 
modes in the harmonic approximation D{uj) for static 
packings of ellipsoidal particles: 1) We compared D(ui) 
obtained from static packings of ellipsoidal particles using 
Perram and Wertheim's exact expression (Eq. [2]) for the 
contact distance between pairs of ellipsoidal particles and 
the Gay-Berne approximation described in Sec. Ill Al 2) 
We calculated D(w) for static packings as a function of 
the tolerance used to terminate energy minimization for 
both the MD and CG methods; and 3) We studied the 
system-size dependence of D(u) in systems ranging from 
N = 30 to 960 particles. 

In Fig. [19l we show that the density of vibrational 
modes D(u)) is nearly the same when we use the Per- 
ram and Wertheim exact expression and the Gay-Berne 
approximation to the contact distance for ellipse-shaped 
particles. D(u>) for static packings of ellipse-shaped par- 
ticles is also not dependent on Vtoij which controls the 
accuracy of the conjugate gradient energy minimization 
(Sec. Ill Bl) . for sufficiently small values. Our calculations 
in Fig. [19] (b) also show that D(u>) is not sensitive to the 
energy minimization procedure (i.e. MD vs. CG) for 
small values of the minimization tolerance i^toi- 

In addition, key features of the density of vibrational 
modes are not strongly dependent on system size. For 
example, in Fig. [20l we show D(oj) for ellipses in 2D at 
aspect ratio a — 1.05 and compression A<fi — 10~ 7 over a 



(20) 



range of system sizes from N - 
D(ui) at fixed system size N 



30 to 960. (For reference, 
240 and A</> = 10~ 8 over 
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FIG. 17: The contribution of the translational degrees of free- 
dom T to each eigenvector m of the dynamical matrix versus 
frequency uj in packings of ellipses in 2D at Acj> — 10 -7 . (a) 
shows data for aspect ratios a = 1.05 (black solid), 1.2 (red 
dashed), 1.5 (green dash-dash-dot), 2.0 (blue dash-dot), and 
4.0 (purple dot-dot-dash) and (b) shows data for aspect ratios 
a = 1.001 (black solid), 1.002 (red dashed), 1.005 (green dot- 
dot-dash), 1.01 (blue dash-dot), 1.02 (purple dot-dot-dash), 
and 1.05 (cyan dotted). The inset to (b) shows (T) averaged 
over modes in the lowest (squares) and intermediate (circles) 
frequency regimes. The solid line has slope 1.5. 




FIG. 18: The contribution of the translational degrees of 
freedom T to each eigenvector m of the dynamical matrix 
versus frequency uj in packings of prolate ellipsoids in 3D at 
A(j> — 10 -6 . (a) shows data for aspect ratios a = 1.01 (black 
solid), 1.05 (red dashed), 1.1 (green dash-dash-dot), 1.2 (blue 
dash-dot), and 1.5 (purple dot-dot-dash) and (b) shows data 
for aspect ratios a — 1.001 (black solid), 1.002 (red dashed), 
1.005 (green dash-dash-dot), 1.01 (blue dash-dot), 1.02 (pur- 
ple dot-dot-dash), and 1.05 (cyan dotted). The inset to (b) 
shows (T) averaged over modes in the lowest (squares) and 
intermediate (circles) frequency regimes. The solid line has 
slope 1.25. 



a range of aspect ratios is shown in Fig. [71) D(oS) in the 
low and intermediate frequency bands and plateau region 
overlap for all system sizes. The only feature of D(oj) 
that changes with system size is that successively lower 
frequency, long wavelength translational modes extend 
from the plateau region as system size increases. In the 
large system-size limit N > (4>— 4>j)~ 2 , which we do not 
reach in these studies, the lowest frequency modes will 



scale as D(lu) ~ uj d 1 . 

IV. CONCLUSIONS 

We performed extensive numerical simulations of static 
packings of frictionless, purely repulsive ellipses in 2D 



14 



3 
2 

1 1 
Q 

rf 

en 
o 

-1 
-2 
-3 





(a) 


: j 




i . i . i . i . 





-4 



-3 -2 
log 10 co 




FIG. 19: (a) Density of vibrational modes in the harmonic 
approximation D{ui) for N — 30 ellipses with a — 1.05 at 
A(f> = 10 -7 using the Perram and Wertheim exact contact dis- 
tance between pairs of ellipses with CG energy minimization 
tolerance Vtoi = 10 -16 (green dot-dashed) and Vtoi = 10~ 8 
(blue dash-dash-dotted) or the Gay-Berne approximation 
with Vtol = 10" 16 (black solid) and Vtol = 10" 8 (red dashed), 
(b) D(u) for N = 12 ellipses with a = 1.05 at A(j> = 10~ 5 
using the Perram and Wertheim exact contact distance with 
CG energy minimization tolerance Vtoi = 10 -16 (solid black), 
and MD energy minimization tolerance K to \ = 10~ 16 (red 
dashed) and 10" 24 ( green dotted). 



and prolate ellipsoids in 3D as a function of aspect ratio 
a and compression from jamming onset A<f). We found 
several important results that highlight the significant 
differences between amorphous packings of spherical and 
ellipsoidal particles near jamming. First, as found pre- 
viously, static packings of ellipsoidal particles generically 
satisfy z < z; so |14l - [l7| ; i.e. they possess fewer contacts 
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FIG. 20: Density of vibrational modes in the harmonic ap- 
proximation D(uj) for ellipses in 2D with aspect ratio a — 
1.05, A<t> = 10" 7 , and system size N = 30 (black solid), 120 
(red dashed), 240 (green dot-dashed), and 960 (blue dash- 
dash-dotted). 



than the minimum required for mechanical stability as 
predicted by counting arguments that assume all contacts 
give rise to linearly independent constraints on particle 
positions. Second, we decomposed the dynamical matrix 
M = H — S into the stiffness H and stress S matrices. 
We found that the stiffness matrix possesses N(zi so — z) 
eigenmodes eo with zero eigenvalues over a wide range 
of compressions (A(f> > 0). Third, we found that the 
modes eo are nearly eigenvectors of the dynamical matrix 
(and the stress matrix — S) with eigenvalues that scale as 
cAcf>, with c > 0, and thus finite compression stabilizes 
packings of ellipsoidal particles [l^]. At jamming onset, 
the harmonic response of packings of ellipsoidal parti- 
cles vanishes, and the total potential energy scales as <5 4 
for perturbations by amplitude S along these 'quartic' 
modes, eo. In addition, we have shown that these results 
are robust; for example, the density of vibrational modes 
D{ui) (in the harmonic approximation) is not sensitive to 
the error tolerance of the energy minimization procedure, 
the system size, and the accuracy of the determination 
of the interparticle contacts over the range of parameters 
employed in the simulations. 

These results raise several fundamental questions for 
static granular packings: 1) Which classes of particle 
shapes give rise to quartic modes?; 2) Is there a more 
general isostatic counting argument that can predict the 
number of quartic modes at jamming onset (for a given 
packing-generation protocol)?; and 3) Are systems with 
quartic modes even more anharmonic [3~ij than packings 
of spherical particles in the presence of thermal and other 
sources of fluctuations? We will address these important 
questions in our future studies. 
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Appendix A: Scaling Behavior of the Total Potential 
Energy 



The scaling behavior of AV/N (shown in Figs. [15] 
and I16|) as a function of the amplitude 6 of the per- 
turbation along the eigenmodes of the dynamical matrix 
is valid only when the original contact network of the 
perturbed static packing does not change. As shown in 
Fig. [5TJ AV/N does not obey the power-law scaling de- 
scribed in Eq. [19] when new interparticle contacts form. 
Note that changes in the contact network are more li kely 
for systems with a ~ 1 as shown previously in Ref. (34j . 
In a future publication, we will measure the critical per- 
turbation amplitude below which new contacts do not 
form and existing contacts do not change for each mode 
k. This work is closely related to determining the non- 
linear vibrational response of packings of ellipsoidal and 
other anisotropic particles. 



Appendix B: Dynamical Matrix for Ellipse-shaped 
Particles 



In this Appendix, we provide explicit expressions for 
the dynamical matrix elements (Eq. lll[) for ellipse-shaped 
particles that interact via the purely repulsive linear 
spring potential (Eq. [1]). The nine dynamical matrix el- 
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and the nine dynamical matrix elements for i = j are 
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where ^ is the polar angle defined in Fig. [SI I = dlfm, 
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FIG. 21: Change in the total potential energy AV/N when 
we (a) do not allow the system to gain contacts or (b) al- 
low the system to gain contacts versus the amplitude of the 
perturbation 5 along several 'quartic' modes (mode 17: black 
solid, mode 23: red dashed, mode 38: green dot-dot-dashed, 
and mode 78: blue dash-dash-dotted) from a static pack- 
ing of N = 240 ellipse-shaped particles at A<j> = 10 -8 and 
a = 1.002. The dotted (dot-dashed) line has slope 4 (2). (c) 
The number of new contacts N' c that differ from the original 
contact network as a function of S for each mode in (a) and 
(b). The vertical lines indicate the S at which the first new 
contact forms for each mode. 
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